Corrections to the neutrinoless double-beta-decay operator in the shell model 
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We use diagrammatic perturbation theory to construct an effective shell-model operator for the 
neutrinoless double-beta decay of 82 Se. The starting point is the same Bonn-C nucleon-nucleon 
interaction that is used to generate the Hamiltonian for recent shell-model calculations of double-beta 
decay. After first summing high-energy ladder diagrams that account for short-range correlations 
and then adding diagrams of low order in the G matrix to account for longer-range correlations, 
we fold the two-body matrix elements of the resulting effective operator with transition densities 
from the recent shell-model calculation to obtain the overall nuclear matrix element that governs 
the decay. Although the high-energy ladder diagrams suppress this matrix element at very short 
distances as expected, they enhance it at distances between one and two fermis, so that their overall 
effect is small. The corrections due to longer-range physics are large, but cancel one another so 
that the fully corrected matrix element is comparable to that produced by the bare operator. This 
cancellation between large and physically distinct low-order terms indicates the importance of a 
reliable nonperturbative calculation. 



I. INTRODUCTION 

Neutrinoless double-beta decay is a very slow lepton- 
numbcr- violating nuclear process that occurs if neutrinos 
are their own antiparticles. An initial nucleus (Z,A), 
with proton number Z and total nucleon number A, de- 
cays to (Z + 2, A), emitting two electrons in the process 

The neutrino masses and mixing matrix figure promi- 
nently in the decay. The rate, assuming that the process 
is mediated by the exchange of a light virtual neutrino, 
is 

Pi/ 2 ] _1 = G 0u (Q, Z) |M „| 2 (m 0f} ) 2 , (1) 

where Q is the energy difference between the initial 
and final nuclei, Z is the charge of the initial nucleus, 
Gq v (Q, Z) is a tabulated phase-space measure, M „ is 
the nuclear matrix element to which we turn shortly, and 
mpff is a- linear combination of neutrino masses: 
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rnw = \J2 m * U ek\ ■ (2) 

k=l 

In this last equation, toj, is the mass of the k th neutrino 
(these mass eigenstates are linear combinations of the 
electron-, mu- and tau-neutrinos) and U e k is the element 
of the unitary mixing matrix that connects that neutrino 
to the electron neutrino. The quantity mpp is what ex- 
perimenters want to extract from the measured decay 
rate. They cannot do so, however, without knowing the 
matrix clement Mq v , which must be calculated in some 
nuclear model. 

Most calculations of M$ v are done cither in the 
neutron-proton Quasiparticle Random Phase Approxi- 
mation (QRPA) or in the shell model. The two methods 
have complementary virtues. The QRPA includes many 



single-particle orbit als outside a relatively small "inert" 
core, but limits itself to a particular kind of correlation. 
The shell model includes arbitrary complicated correla- 
tions, but only among a few single-particle orbitals out- 
side a larger inert core. The current predictions of the 
two models, after a recent shaking out period, show the 
QRPA matrix clement exceeding that of the shell model 
by factors of up to about two in the lighter isotopes such 
as 76 Ge and 82 Se, and somewhat less in the heavier iso- 
topes 

Which kind of calculation is closer to the truth? Are 
there important effects that escapes both models? To 
find out, one has to correct one or both to account for 
omitted physics. Although it is possible to add miss- 
ing correlations to the QRPA, it is not easy to do so 
systematically because several different uncontrolled ap- 
proximations — BCS pairing, a phcnomcnological inter- 
action, the quasiboson approximation, etc. — are part 
of the method. By contrast, because the shell model in- 
cludes all correlations within a well-defined subspace of 
the full Hilbert space (the space generated by valence 
particles occupying a few single-particle states), there is 
a systematic procedure for adding the effects of states 
outside that space [HE]. 

While the procedure as usually implemented is pertur- 
bative in a renormalized residual nuclear interaction (the 
G matrix) and not always reliable for that reason, it of- 
ten works well enough, particularly if followed by some 
modest adjustment to data. Practitioners have long used 
such an approach to obtain good effective interactions [5] , 
but have never applied the same techniques to obtain an 
effective double-beta-decay operator. Instead, they typ- 
ically modify the bare operator phenomcnologically, e.g. 
through the reduction of the axial-vector coupling con- 
stant gA (suggested by studies of single-beta decay) or the 
use of a prescription [7] to treat the short-range nucleon- 
nucleon repulsion that is not present in shell-model wave 
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functions 1 . Not surprisingly, it is difficult to assess the 
reliability of such approximations. 

In this paper, therefore, we apply the same techniques 
used to construct effective shell- model interactions to the 
decay operator itself. Section II below contains a brief 
description of the matrix element we try to calculate, 
and Section III a description of our procedure for renor- 
malizing the corresponding operator. In Section IV we 
present the results of our calculation, which combines the 
renormalized operator with shell-model transition densi- 
ties from the authors of Ref. [TO] for the decay of 82 Se 
to 82 Kr. (Densities for other decays, e.g. of 76 Ge, are 
not currently available.) Although we cannot be fully 
confident in our perturbative result, our initial steps — 
accounting for short-range correlations through the gen- 
eration of a G matrix and an analogous corrected decay 
operator — are nonperturbative and trustworthy. And 
even if the low-order perturbation theory we employ sub- 
sequently is not accurate, it should tell use whether we 
might expect significant renormalization in a fully non- 
perturbative treatment. 



II. FORM OF BARE M 0v 

A precise expression for the matrix element is com- 
plicated, but with a few approximations that induce an 
error of with an error of less than 30% |H[5], we can write 
M 0v as 



GT 9V M F 



(3) 



with gy and gA the vector and axial-vector coupling con- 
stants, and 



Ml = (f\J2H(r ab ,E) T +T+\ 



M, 



GT 



(f\J2H(r ab ,E)a a -a b T+Tj 



(4) 



(5) 



Here \i) and |/) are the initial and final nuclear ground 
states, a, b label nucleons, E is an average excitation en- 
ergy, and H is a "neutrino potential" , given by 



H(r,E) = 



2R 



dq 



sin qr 



q + E-{Ei+Ef)/2 



(6) 



The quantity R is the nuclear radius, inserted to make 
the matrix element dimcnsionless. Since our work is ex- 
ploratory, we use the relatively simple forms in Eqs. Q 
and ([5| in most of what follows, though we also discuss 
corrections due to nucleon form factors and forbidden 
terms in the weak nuclear current. 



1 The use of the Unitary Correlation Operator Method, e.g, in 
Ref. [HI [9| , is more than a prescription, but the method has not 
yet been consistently applied to both the decay operator and the 
interaction. 



III. CONSTRUCTING AN EFFECTIVE 
SHELL-MODEL DECAY OPERATOR 

A. Formalism and diagrams for two-body operators 

Diagrammatic effective-operator theory has a long his- 
tory in nuclear physics. References [IT] and j4] are early 
reviews and Ref. is a more recent one. The theory is 
based on the division of the many-body Hilbert space into 
the shell-model space P of particles occupying several de- 
generate or quasidegenerate orbitals (usually eigenstates 
a harmonic-oscillator potential U), and the rest of the 
Hilbert space Q. One begins by defining operators with 
the same names that project onto these spaces: 



i eP 



other i 



ith 



P z = P . 



Q 2 = Q, 



PQ = QP = 0. 



(8) 



Next one defines an effective Hamiltonian that when act- 
ing on the P-space projection of an eigenstate l^o) gives 
back that projection with the correct eigenvalue: 



H oS (E a )P\V a ) = E a P\^ a ) 



(9) 



Similarly, for any "bare" operator M. (where "bare" 
means "acting in the full model space"), we can define 
an effective operator M e s that acts only in the P space, 
with matrix elements related to those of the bare opera- 
tor by: 



(y a \PM cH P\V b ) 



(10) 



It is straightforward to show that the effective Hamil- 
tonian can be represented in energy-dependent form as a 
solution to the Bloch-Horowitz equation [TJ] 



H cS (E) = PHP + PHQ 



1 



E-QH 



QHP, (11) 



with the full wave function a solution to the associated 
equation 

1 



\HE))=Z 1 



E-QH 



QH PME)) , (12) 



Here E is the energy of the eigenstate and Z is a normal- 
ization factor. One can remove the explicit dependence 
on energy by treating the residual interaction V — U , 
which couples the P and Q spaces, as a perturbation. 
When Eqs. (Ill and (12 1 are solved order by order, the 



result is a series of valence- linked Goldstone diagrams for 
the matrix elements matrix elements of iJ ff and M. e a 

The diagrams are very much like those for the binding 
energy, but have open lines at each end to represent the 
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valence single-particle states on the right and left sides 
of the effective operators. The use of this representation 
in a linkcd-cluster expansion forces the introduction of 
"folded" diagrams, in which intermediate states have zero 
excitation energy. The mostly low-order diagrams we 
consider here, however, will not have folds. 

Some recent work on effective interactions |13j has used 
Mow k as a starting point, with high-energy states effec- 
tively integrated out at the beginning. But because we 
want to calculate the contributions of such states to the 
effective decay operator, we need a method that treats 
them explicitly. Our starting point, therefore, is the same 
as in the traditional treatment of nuclear matter: we de- 
fine a nonperturbative G matrix as the sum of the two- 
particle ladder diagrams displayed in Fig. [T] below. The 
ruled lines in the figure indicate high-energy states, lying 
well above the shell-model single-particle space (we will 
vary the exact amount by which they are above). The G 
matrix is thus defined not only between two-body states 
in the valence space but also between states in a larger 
model space that contains several higher shells. (See Ref. 
5 for details on this "double partitioning" of the Hilbert 
space.) 



G 



+ » 




- + 



+ 



FIG. 1: Diagrammatic construction of a G matrix, the first 
step in an effective interaction. Ruled lines correspond to 
high-energy particle states in the doubly partitioned space. 

Although there are an infinite number of ladder dia- 
grams, the sum can be carried out indirectly, e.g. through 
the solution of the Bethe-Goldstone equation. The fa- 
miliar idea underlying the infinite sum is that the hard 
short-range core that makes the nucleon-nuclcon inter- 
action intractable can be treated exactly, at least at the 
two-body level, by nonperturbatively admixing into the 
wave function intermediate two-particle states with arbi- 
trarily high-energy. The effective low-energy interaction 
G[V] that results has a soft core because the effects of 
short-range physics have already been accounted for in 
the ladder sum. The argument V in the G matrix is 
meant to reinforce the fact that G depends on the "bare 
interaction" V. 

After this nonperturbative construction of the G ma- 
trix, one can use perturbation theory in G to add the 
effects of states that are at low energy but still outside 
the valence space, i.e. in the intermediate space of the 



'\/V\/V 






FIG. 2: Low-order diagrams (in G) for the effective interac- 
tion H c ft . 



double-partitioned set. The diagrams in Fig. [2] include 
all such effects up to second order in G (with the excep- 
tion of those produced by tadpole and one-body graphs, 
which are commonly omitted), and some third-order ef- 
fects. In this figure, the upward-going lines represent 
low-lying particle states (including the valence levels, as 
long as they do not lead to intermediate denominators 
with zero energy) in the gap between the Fermi surface 
and the high-energy levels . Downward going lines rep- 
resent "hole" states that correspond to the vacating of 
levels below the shell-model space. 

Typically, more complicated graphs, including the 
folded ones, are included alongside the graphs in the fig- 
ure. Even then, problems with convergence, three- and 
higher-body operators that are generally too complicated 
to include, etc., mean that the resulting interaction often 
must be modified phenomenologically, especially in the 
monopolc-monopole channel [Hj Sometimes, how- 

ever, the perturbation theory by itself is enough to pro- 
duce a pretty good interaction [5J |T3] . The recent shell- 
model calculations of double-beta decay in Ref. jTU] were 
based on a tuned version of such an interaction. 

A similar procedure can be followed to evaluate the 
matrix elements of any two-body operator 2 Jv[ . The ex- 



pansion of equation (12 1 leads to a set of diagrams for 
the effective operator in Eq. (10 1 in which a horizontal 



line representing the bare operator appears once along- 
side an arbitrary number of interaction lines |U [TTJ . The 
denominator of Eq. (10 1 gives rise, in addition, to norm 



and overlap diagrams that determine a special basis in 
which the effective operator should be represented. We 
have, however, evaluated the most important of these di- 
agrams (they are given by the derivative with respect to 
the unperturbed energy of the effective-interaction dia- 
grams) and found them to contribute at most a few per- 
cent. We therefore won't include their complicated but 
small effects here. 



2 Previous work has focused almost entirely on one-body opera- 
tors, however. 
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To obtain the effective decay operator, we begin by 
summing all diagrams with two particles excited to high 
energies in which one horizontal line in each diagram is 
the operator M. rather than the interaction. We denote 
the result of this nonperturbative sum, which is com- 
pletely analogous to the G matrix, by .Mhigh- Sequences 
of 1, 2, . . .interaction lines either before or after the op- 
erator insertion can be separately summed, i.e. replaced 
by G matrices (or, more, precisely, the similar ladder sum 
G for which the outgoing states are high- lying). If M. is 
one of the double-beta operators appearing in equations 
|3]) , then two neutron lines become proton lines whenever 
it acts, and the ladder sum reduces to the four diagrams 
in Fig. [3j Thus, the solid lines (red online) in the fig- 
ure represent neutrons and the dotted lines (blue online) 
represent protons. 

Since these diagrams involve only T = 1 states, the 
their sum (with small Coulomb effects neglected) can be 
calculated simply from a G-matrix code through the trick 



■Mhigh = -^G[V T =i 



XM'\ 



(13) 



where M! is the charge-conserving version of the charge- 
changing operator in Eqs. Q or ([5]), obtained by remov- 
ing the isospin raising operators. In other words, we 
can calculate matrix elements of .Mhigh by computing 
the G matrix corresponding to the interaction V + \M! '. 
The derivative filters out all graphs except for those that 
have a single double-beta line replacing an interaction 
line. The difference between the matrix elements of M. 
and .Mhigh give us a rigorous measure, at least for two- 
valence-nucleon systems, of the effects of short-range cor- 
relations in double-beta decay (up to the few percent due 
to norm diagrams). 

Having constructed .Mhigh to include short-range two- 
body correlations, one can use it together with the G ma- 
trix to calculate the additional renormalization from low- 
lying excitations. We do so by replacing one G-matrix 
line in each of the diagrams in Fig. [2] by .Mhigh- The 
resulting diagrams for M c fr, all first order in G except 
for second-order ladders, appear in Fig. [4j Short-range 
correlations are included at every vertex through the use 
of G and Mhigh in place of V and M. Most of these 
diagrams can be calculated through the trick in equation 
(131. Only the core-polarization graphs (the last two 



in the figure) must be treated explicitly. Those graphs 
are essentially different from the corresponding effective- 
interaction graph in the identities (neutron or proton) of 
the particles and holes involved. One cannot see the dif- 
ference explicitly in our figures because Fig. [2] does not 
distinguish between neutrons and protons. 



B. Using the operator in 82 Se 



Mush 



+ 




+ 



+ 




FIG. 3: (color online) An effective double-beta operator that 
accounts for short-range correlations in the nuclear wave func- 
tion. Solid lines (red online) are neutrons and dotted line 
(blue online) are protons. The symbol G represents the ex- 
tension of G to the Q space. 



more diagrams). The Bonn-C nucleon-nucleon interac- 
tion |15j was the starting point for the ladder-diagrams 
that enter the G matrix. After summing all diagrams 
through third order in G and the folded diagrams based 
on those, the authors adjusted the interaction by fitting 
certain components to spectra. 

To be as consistent as possible with the calculations 
of Ref. [TO] in constructing our effective decay operator 
for 82 Se, we use the same Bonn-C interaction, the same 
valence space (/5/2PS9/2)! the same oscillator parameter 
(b = 82 1 / 6 MeV), and the same average energy (E — 
10.08 MeV) as that reference throughout. After calcu- 
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The shell-model calculations of Ref. [lOj used an inter- 
action that was constructed largely through the effective- 
interaction theory discussed above (though it included 



FIG. 4: (color online) Low order diagrams (in G) contributing 
to the final effective double-beta operator. The thick dashed 
line is .Mhigh, expressed diagrammatically in Fig. [2] 
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lating the G matrix and its extension to the Q space, we 
evaluate .Mhigh in Fig. [3] and the more involved effective- 
operator diagrams in Fig. |4j Our only modification to the 
standard procedure is to prohibit intermediate particles 
in diagrams (i) and (j) from occupying levels that are 
essentially full in 82 Se and 82 Kr; these Pauli-forbidden 
contributions, if included, would be canceled by higher- 
order diagrams that we do not evaluate here. Finally, we 
combine the matrix elements of the effective two-body 
operator (at several stages of approximation) with the 
65 independent two-body ground-state-to-ground-state 
transition densities from the shell-model calculation of 
Ref. |10j to obtain a transition matrix element for the 
decay 82 Se — > 82 Kr. 



IV. RESULTS 

A. High-energy states and short-range correlations 

Before looking at long-range corrections, we report the 
effects of the ladder diagrams shown in Fig. [3] The treat- 
ment of short-range correlations these diagrams repre- 
sent is completely well defined; one knows exactly what 
it includes and what it omits, and there is no double 
counting. To look at the spatial structure of the correla- 
tions we define a two-body double-beta correlation func- 
tion C _ (r) by making the substitution H(r a b,E) — > 
H(r a b,E) S(r — r b) in the Gamow- Teller transition op- 
erator. 

Figure [5] displays the results for the decay of 82 Se 
with the shell-model transition densities mentioned 
above, when the boundary between "high-energy" single- 
particle states (ruled in the diagrams) and lower-energy 
states lies 4 hcu above the valence fp shell. The results 
change only very slowly as the boundary is moved up 
from that point. The function in the top panel, labeled 
with the subscript 0, corresponds to the simple GT opera- 
tor of Eq. (|5| ; the function in the bottom panel includes 
modifications to that operator from the weak nucleon 
form factors and higher-order terms in the weak current 
(see, e.g., Refs. [3] and [16] for definitions). 

Both panels show the suppression of short-range con- 
tributions by the ladders, though the suppression is 
weaker in the lower panel because the nucleon form fac- 
tors cut out some short-distance neutrino exchange by 
themselves. But in neither case is the matrix element 
reduced very much by the correlations: for the simple 
operator the reduction is about 8%, and for the full op- 
erator it is less than 3%. The reason, as the figure shows, 
is that probability density is shifted from very short dis- 
tances to around r = 1 fm, where the function H(r, E) 
is still large. 

This behavior is similar to that found in Ref. [To] , 
which constructs correlated wave functions by solving 
coupled-cluster equations at the one and two-body level. 
There are some differences between that treatment of 
short-range effects and ours: the coupled-clusters method 




FIG. 5: (color online) The radial distribution C (r) that, 
when integrated, produces the Gamow- Teller 0v matrix ele- 
ment for the decay of 82 Se. The top panel corresponds to 
the simple operator given in Eq. (JsJ and the bottom panel 
includes momentum-transfer dependent form factors and for- 
bidden operators (see text). The solid lines correspond to the 
bare operator and the dashed lines to the effective operator 
generated by summing the high-energy ladders in Fig. [4] 



includes Hartree-Fock-like effects, sums hole-hole ladders 
as well as particle-particle ladders, and has no double 
partitioning. Nevertheless, both methods include much 
of the same physics and should yield similar results. To 
test the similarity we repeated our calculation with the 
Argonne V18 interaction used if Ref. [TB] in place of 
the Bonn-C interaction. Figure [6] shows the results for 
the simple-operator GT correlation function discussed 
above along with the corresponding coupled cluster re- 




FIG. 6: (color online) Same as the top panel of Fig. [5] but 
with the dashed line from a calculation with the Argonne 
V18 potential in place of Bonn-C. The dotted line, nearly 
indistinguishable from the dashed line, comes from a fit to 
coupled-cluster results [TS] with the Argonne potential. 
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TABLE I: Renormalization of Mq^ (from equation |5|) for 
the decay 82 Se — > 82 Kr. 



boundary 


bare 


(a) 


(a)-(f) 


(a)-(h) 


all 


4Hlo 


3.33 


3.07 


4.15 


5.38 


3.05 


5hcu 


3.33 


3.06 


4.17 


5.39 


3.15 


6hui 


3.33 


3.05 


4.16 


5.39 


3.21 


7huj 


3.33 


3.06 


4.17 


5.39 


3.28 


8huj 


3.33 


3.06 


4.17 


5.39 


3.35 



suits, which we generated by using the phenomenological 
J astrow- function fit reported in Rcf. [IB] . The two cor- 
relation functions are almost indistinguishable. These 
results, alongside those of Ref. 0E], mean that the phe- 
nomenological J astrow function from Ref. [7], used re- 
flexively for a long time, almost certainly overestimates 
the quenching due to short-range correlations. 

B. Longer-range correlations 

When the correlations induced by the diagrams in Fig. 
|3]are added to the short-range effects discussed above, the 
matrix element changes further. Results for the Gamow- 
Teller part of the matrix element (equation ^ , without 
form factors or forbidden currents) appear in table [I] 

The table successively adds the results from three 
classes of diagrams: the high-energy ladders discussed 
above (labeled (a) in Fig.[4|, the diagrams with the high- 
energy ladders embedded in lower energy ladders (labeled 
(b) - (f)), the 4-particle 2-hole diagrams (g) and (h), and 
finally the core-polarization diagrams (i) and (j). Each 
row corresponds to a different boundary (measured in 
hui from the fp shell) between the "high-energy" parti- 
cle levels, denoted by ruled lines in the figures, and the 
lower-energy levels. All the diagrams except the last two 
in Fig. [4] are insensitive to this boundary if it is above 
about 4 huj. The core polarization graphs, by contrast 
still have not converged at 8 (or even 9) Tmj. We are 
unable to carry the calculation beyond that point. Core- 
polarization graphs in the effective interaction are noto- 
rious for converging very slowly, sometimes taking 20 or 
more fuv |17j . 

The table shows several things. First, the short-range 
correlations, as discussed previously, damp the bare ma- 
trix element by about 8%. Second, the ladder and 4- 
particle- 2-hole contributions then increase the matrix ele- 
ment by about 75%. These graphs contain pairing matrix 
elements that promote particles into unoccupied levels, 
unblocking the transition. Finally, the core-polarization 
diagrams decrease the matrix element so that after sum- 
ming particle-hole configurations up to 8 hto our matrix 
element is only marginally bigger than the bare version. 



The core-polarization graphs contain the neutron-proton 
interaction, the correlations from which are known to 
counteract the effects of pairing. Cutting off the sum 
at 8 hut, however, probably exaggerates the size of this 
counteraction; the bulk of the effect comes from low en- 
ergy levels, and contributions from higher-energy levels 
actually increase the matrix element again. Thus, our 
full result (in the column labeled "all") grows with the 
boundary between low and high-energy states, and is still 
growing at our maximum value. 

All these statements remain true when true when we 
include the Fermi term in Eq. Q — resulting in a total 
M 0v of 3.95 at 8 huj vs. the bare value 3.78 — or add the 
effects of form factors. We have not included forbidden 
currents in the full calculations, but do not expect them 
to change the pattern reported in the table. 

As just noted, the effects of the pairing and neutron- 
proton correlations cancel each other to a significant ex- 
tent. The cancellation resembles what happens in the 
QRPA, which includes a portion of the effects calculated 
here through the use of a relatively large single-particle 
space. When all is said and done, our result, is a bit 
larger (and continuing to grow at 8hu) than result at the 
bottom of column (a). That number is the bare matrix 
element corrected for short-range correlations, like the 
results reported in Refs. [2] and [10]. The use of our effec- 
tive operator thus improves the agreement between the 
shell model and QRPA, though the convergence issues 
keep us from saying by exactly how much. One would 
expect something similar for 76 Ge, which has only six 
fewer nucleons and is treated in the same model space. 

Of course, the large corrections from individual graphs 
mean that higher order corrections may be sizeable as 
well. One cannot, therefore, take the results of this first- 
order estimate too seriously. In particular, it would not 
be correct to conclude that the near cancellation between 
pairing and neutron-proton correlations persists to higher 
order and/or when many-body diagrams are included. 
The large individual components indicate that the con- 
tributions of states outside the model space are signifi- 
cant and may or may not cancel one another in a better 
calculation. A more reliable nonperturbativc evaluation 
of the corrections is very important, and we are working 
in that direction. 
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